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ABSTRACT 

We report the detection of high-energy 7-ray emission from two starburst 
galaxies using data obtained with the Large Area Telescope on board the Fermi 
Gamma-ray Space Telescope. Steady point-like emission above 200 MeV has 
been detected at significance levels of 6.8<r and 4.8a respectively, from sources 
positionally coincident with locations of the starburst galaxies M82 and NGC 253. 
The total fluxes of the sources are consistent with 7-ray emission originating from 
the interaction of cosmic rays with local interstellar gas and radiation fields and 
constitute evidence for a link between massive star formation and 7-ray emission 
in star-forming galaxies. 

Subject headings: Galaxies: individual (M82, NGC 253) — gamma rays: obser- 
vations — cosmic rays — radiation mechanisms: non-thermal 

1. Introduction 

Cosmic rays are believed to be accelerated by supernova remnant shocks that are formed 
when a star explodes (Ginzburg & Syravotskii 1964; Hayakawa 1969). Observations of 7 rays 
from supernova remnants in the Milky Way would apparently offer the best opportunity to 
identify the sources of cosmic rays, but cosmic-ray diffusion throughout the Galaxy results in 
a bright 7-ray glow, making it difficult to attribute 7 rays to cosmic-ray electrons, protons or 
ions accelerated by Galactic supernova remnants. Direct evidence for the sources of cosmic 
rays is therefore still lacking. 

The supernova remnant paradigm for cosmic-ray origin can also be tested by measuring 
the 7-ray emission from star-forming galaxies. Starburst galaxies, in particular, should have 
larger 7-ray intensities compared to the Milky Way due to their increased star-formation rates 
and greater amounts of gas and dust that reprocess light into the IR, and, with photons, serve 
as targets for 7-ray production by cosmic ray electrons and ions. If the 7-ray production rate 
is sufficiently increased, star-forming galaxies will be detectable by the current generation of 
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instruments, as early estimates (e.g., Volk et al. 1989, 1996; Akyiiz et al. 1991; Paglione et 
al. 1996) and recent detailed models (e.g., Domingo-Santamaria & Torres 2005; Persic et al. 
2008; de Cea del Pozo et al. 2009; Rephaeli et al. 2009; Lacki et al. 2009) predict. 

Here we report the detection of the starburst galaxies M82 and NGC 253 in high-energy 
7 rays from observations with the Large Area Telescope (LAT) on board the Fermi Gamma- 
ray Space Telescope. A description of the analysis of the observations is given in Section 
2. In Section 3, the measured spectra and fluxes are compared with predictions based on 
theories of cosmic-ray origin from supernovae in star-forming galaxies. 

2. Observations and Analysis 

The LAT is a pair-conversion telescope with a precision tracker and calorimeter, a seg- 
mented anti-coincidence detector (ACD) which covers the tracker array, and a programmable 
trigger and data acquisition system. Incoming 7 rays convert into electron-positron pairs 
while traversing the LAT. The directions of primary 7 rays are reconstructed using informa- 
tion provided by the tracker subsystem while the energies are measured via the calorimeter 
subsystem. The ACD subsystem vetoes the great majority of cosmic rays that trigger the 
LAT. The energy range of the LAT spans from 20 MeV to > 300 GeV with an angular 
resolution of approximately 5.1° at 100 MeV and narrowing to about 0.14° at 10 GevQ. Full 
details of the instrument, onboard and ground data processing, and other mission-oriented 
support are given in Atwood et al. (2009). 

The LAT normally operates in a scanning mode (the 'sky survey' mode) that covers 
the whole sky every two orbits (i.e., ~ 3 hrs). We use data taken in this mode from the 
commencement of scientific operations in early-August 2008 to early- July 2009. The data 
were prepared using the LAT Science Tools package. Only events satisfying the standard low- 
background event selection (the so-called 'Diffuse' class events corresponding to the P6V3 
instrument response functions described in Rando (2009)) and coming from zenith angles 
< 105° (to greatly reduce the contribution by Earth albedo 7 rays) were used in the present 
analysis. To further reduce the effect of Earth albedo backgrounds, time intervals when the 
Earth was appreciably in the field of view (specifically, when the center of the field of view 
was more than 43° from the zenith) were also excluded from the analysis. 



1 Angular resolution is denned here as the 68% containment radius of the LAT point spread function 
averaged over the intrument acceptance and including photons which convert in either the thick or thin 
layers of the tracker array. 
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We use all 7 rays with energy > 200 MeV within a 10° radius region of interest (ROI) of 
the optical locations for the galaxies M82 and NGC 253. Detection significance maps for each 
ROI are shown in Fig. [TJ The background model for each ROI includes all LAT-detected 
sources along with components describing the diffuse Galactic and isotropic 7-ray emissions. 
Each map shows a bright and isolated 7-ray excess above the background that is consistent 
with the location of the nominal (optical) position of the respective starburst galaxy. 

The data were analyzed using the LAT Science Tools package (v9rl5p2), which is avail- 
able from the Fermi Science Support Center, using P6V3 post-launch instrument response 
functions (IRFs). These IRFs take into account event pile- up and accidental coincidence ef- 
fects in the detector subsystems that were not considered in the definition of the pre-launch 
IRFs. We used a maximum likelihood fitting procedure (gtlike) to determine the positions of 
the 7-ray sources associated with M82 and NGC 253 (see Tabled]). The angular separation 
between the best-fit location and the core of each galaxy is 0.05° for M82 and 0.12° for 
NGC 253. Systematic uncertainties in the positions due to inaccuracies in the point spread 
function and telescope alignment are estimated to be less than 0.01°. 

We tested the possibility that the sources are spatially extended by fitting two-dimensional 
Gaussian-shaped intensity profiles. The widths and locations of the profiles were adjusted 
and refit over the region in an iterative procedure but we found no significant evidence for 
source extension in our data. We verified these results using a likelihood fitting procedure 
capable of modeling spatially extended 7-ray sources (sourcelike). A comparison between 
the point and extended source hypotheses using this method produces negligible changes 
in detection significance. From our analysis we set upper limits on the angular sizes of the 
emitting regions as 0.18° for M82 and 0.30° for NGC 253 at the 95% confidence level assum- 
ing a two-dimensional Guassian spatial model parameterized by the 68% surface intensity 
containment radius. By comparison, the angular sizes of the galaxies are 0.19° x 0.07° for 
M82 and 0.46° x 0.11° for NGC 253 as measured in the ultraviolet band (Gil de Paz et al. 
2007). The starburst cores of M82 (Volk et al. 1996) and NGC 253 (Ulvestad 2000) have 
an angular extent < 0.01° and cannot be resolved by the LAT. 

Spectral analysis is a separate maximum likelihood calculation for which we have adopted 
the point-source hypothesis and best-fit position determined during the localization and ex- 
tension fitting step. 

Diffuse 7-ray emission from the Milky Way is treated with the Galactic diffuse model 
described within the gllJem_v02.fit file suitable for analysis with the Science Tool gtlike. 
In addition to the spatially structured Galactic diffuse emission, Fermi also observes an 
isotropic diffuse component which includes both extragalactic diffuse 7-ray emission and 
instrumental background from charged particles triggering the LAT. The isotropic diffuse 
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emission has been treated with isotropic_iem_v02.txt. Excepting the sources associated with 
M82 and NGC 253, all individual objects detected by Fermi after 11 months of scientific 
operations within a 10° radius of the best-fit position of each galaxy are also included into 
the background description of each region as distinct point sources. 

We considered alternative associations for the two LAT sources of interest aside from 
M82 and NGC 253 in the CRATES catalog of flat-spectrum radio sources (14467 entries, 
Healey et al. 2007) and the Candidate Gamma-Ray Blazar Survey catalog, CGRaBS (1625 
entries, Healey et al. 2008). Both of these catalogs show high correlation with 7-ray bright 
blazars based on multiwavelength observations. However, there are no likely CRATES or 
CGRaBs objects within the positional uncertainty of either LAT source. Near NGC 253, 
the only source of possible concern is a ~ 40 mJy NVSS (Condon et al. 1998) radio source 
at 1.4 GHz with unknown spectrum. Such a source would be unusually weak by comparison 
with the radio fluxes of LAT blazars. 

We searched for flux variability for each 7-ray source by creating monthly flux histories of 
the total photon flux > 400 MeV arriving from within a circular regions 1° in radius centered 
on the Fermz-determined locations. No flaring events are observed and the x 2 goodness-of-fit 
test is consistent with constant flux for each source (reduced \ 2 = 0.80 and 1.03 for M82 and 
NGC 253, respectively, each with 9 degrees of freedom). Lack of variability is in accord with 
the cosmic-ray origin hypothesis where most of the emission derives from diffuse cosmic-ray 
interactions, though mild variability of 7 rays and radio emission (Kronberg et al. 2000, 
Brunthaler et al. 2009b) might still occur if M82 or NGC 253 had a recent supernova. Large 
amplitude 7-ray variability on short timescales would rule out a cosmic-ray origin of the 7 
radiation. 

Table [T] summarizes the results of the analyses of M82 and NGC 253. The overall 
detection significance is 6.8<r for M82 and 4.8<r for NGC 253. Note that the significance level 
for these moderately hard spectrum sources is based on the number of high energy photons 
compared to the expected background, whereas the flux uncertainty is based on the number 
of such photons, which is not large, and systematic effects. The integral photon fluxes over 
100 MeV are calculated by extrapolation of the fitted spectral models. 

3. Interpretation 

With the nearest luminous starburst galaxies, M82 and NGC 253, detected by the Fermi 
Gamma-ray Space Telescope, we can test long-standing predictions based on the cosmic-ray 
paradigm that diffuse 7-ray emission from star-forming galaxies is produced via cosmic-ray 
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interactions. The distance to M82 is 3.63 ± 0.34 Mpc (Freedman et al. 1994), and distance 
estimates to NGC 253 range from 2.5 Mpc (Turner & Ho 1985, Mauersberger et al. 1996) 
to 3.9 ± 0.37 Mpc (Karachentsev et al. 2003). Vigorous star formation is observed within 
the central several hundred parsecs of these galaxies. Estimates of the supernova explosion 
(SN) rate vary from « 0.08 - 0.3 yr" 1 in M82, to « 0.1 - 0.3 yr" 1 in NGC 253, compared to 
the supernova rate of 0.02 yr" 1 in the Milky Way. Recent studies of M82 find 7 x 1O 8 M 
in atomic H I gas and 1.8 x 10 9 M Q in H 2 gas (Casasola et al. 2004). The central region of 
NGC 253 contains a bar of molecular gas with an estimated mass of 4.8 x 10 8 M Q (Canzian 
et al. 1988), and its total gas content is ~ 60% of the Milky Way's (Boomsma et al. 2005; 
Houghton et al. 1997; Brunthaler et al. 2009b), reflecting active star formation taking place 
in these relatively small galaxies. 

Table [2] gives adopted values of distance d, supernova rate Rsn, total gas mass M Gas , 
7-ray flux F(> 100 MeV), and 7-ray luminosities for M82 and NGC 253, alongside those 
of the Large Magellanic Cloud (LMC) and the Milky Way. The 100 MeV to 5 GeV 7-ray 
luminosity of M82 and NGC 253 is ~ 10 40 erg s _1 , compared to ~ 3 x 10 39 erg s" 1 for the 
Milky Way, and w 4.1 x 10 38 erg s" 1 for the LMC. These galaxies lack active central nuclei 
and so require a different origin for their 7-ray fluxes than from galaxies with supermassive 
black-hole jets. The 7 rays from our Galaxy and the LMC arise predominantly from cosmic 
rays interacting with interstellar gas and radiation fields. The starburst galaxies M82 and 
NGC 253, though having less gas than the Milky Way, have a factor 2-4 greater 7-ray 
luminosity, suggesting a connection between active star formation and enhanced cosmic-ray 
energy densities in star-forming galaxies. 

We examine several possible correlations between total gas mass, supernova rate, and 
7-ray luminosity of these four galaxies as illustrated in Fig. [2] (cf. Pavlidou & Fields 2001, 
for local group galaxies). In the left-hand panel, we find a poor correlation between 7- 
ray luminosity and gas mass, and a weak linear correlation between 7-ray luminosity and 
supernova rate. Models that attribute the 7 rays to cosmic-ray processes depend both on 
enhanced cosmic-ray intensities, which depends on the supernova rate, and large quantities 
of target gas, suggesting that the 7-ray luminosity is proportional to the product of the 
total supernova rate and gas mass, as shown in the right-hand panel of Fig. [2J Note that 
while the detection of galaxies in this sample is flux-limited, the measured gas masses and 
supernova rates for all galaxies are not, so that the dependence of 7-ray luminosity on these 
parameters reflect underlying physical relationships rather than sensitivity effects. Although 
the sample size is small, this result argues in favor of a scaling of 7-ray luminosity according 
to expectations from the hypothesis that the emission is produced by cosmic-ray interactions. 

Evaluation of the dependence of 7-ray luminosity on galaxy properties is complicated, 
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however, by star formation rates that depend on location in the galaxy. Radio and infrared 
observations reveal that the starburst activity in M82 and NGC 253 takes place in a relatively 
small central region, radius ~200 pc for both M82 (Volk et al. 1996) and NGC 253 (Ulvestad 
2000), so that the distribution of the cosmic rays in the galaxies is probably not uniform. In 
cases where 7-ray emission can be resolved, as for the Milky Way, this can be seen directly 
(Dragicevich et al. 1999). For instance, 7-ray emission from the LMC is mostly produced in 
the star-forming region 30 Doradus, and does not simply trace star formation and total gas 
mass (Abdo et al. 2009a). 

Theoretical predictions, despite using different assumptions and treating the processes 
with varying levels of detail, are largely consistent with the detected integral flux of M82 
(e.g. Volk et al. 1989; Akiiz et al. 1991; Persic et al. 2008; de Cea et al. 2009) and NGC 
253 (e.g. Paglione et al. 1996; Domingo-Santamaria & Torres 2005; Persic et al. 2008). Fig. 
[3] shows the predicted and observed spectra. In the case of NGC 253, the predicted photon 
flux (> 100 MeV) is 2.3 x 10~ 8 photons cm -2 s" 1 (Domingo-Santamaria & Torres 2005) and 
2 x 10~ 8 photons cm -2 s" 1 (Persic et al. 2008). For M82, the predicted photon flux (> 100 
MeV) is between 2.6 x 10~ 8 and 8.3 x 10^ 9 photons cm~ 2 s^ 1 (de Cea et al. 2009) due to 
systematic uncertainties in model parameters; and « 10~ 8 photons cm -2 s~ x (Persic et al. 
2008). Furthermore, extrapolation of the best-fit power-law spectral model at GeV energies 
provides a smooth connection to flux densities of M82 reported at TeV energies (Acciari et 
al. 2009). Although not highly constraining due to the faintness of M82 in the GeV band, 
the fitted spectrum suggests that a single physical emission mechanism dominates from GeV 
to TeV energies. The relationship between the GeV and TeV emission for NGC 253 is less 
clear given the current data. Also, note that the inner starburst region of NGC 253 has 
about a factor of 3 less radio flux than that of M82 at 1.4 GHz, consistent with the galaxy 
being less luminous in 7 rays (M82, Klein et al. 1988; NGC 253, Carilli 1996). 

The star-forming galaxy contribution to the extragalactic 7-ray background (EGB) can 
be estimated by writing the EGB intensity as el** = RuCpbL^/4n, where the Hubble radius 
.Rh — 4200 Mpc for a Hubble constant of 71 km s _1 Mpc -1 , £ ~ 3 - 10 is a cosmological 
factor accounting for more active star formation at redshift z > 1, and p = ^3/(1000 Mpc 3 ) 
is the local space density of normal and star-forming galaxies. The factor b = 0.4 corrects 
for the intensity at 100 MeV given the > 100 MeV luminosity. Writing L 7 = 10 40 L 40 erg s _1 
gives el s j = 3.5 x 10~ 10 b(psL i0 erg cm~ 2 s -1 sr -1 . For L* galaxies like the Milky Way, p 3 = 3 
- 10, and for starburst galaxies like M82 and NGC 253, p^ is an order of magnitude smaller 
(e.g., Scoville 1992). At 100 MeV, a diffuse intensity of e/f GB (100 MeV) ^ 2.4 x 10 -9 erg 
cm -2 s" 1 sr -1 was measured with EGRET (Sreekumar et al. 1998), similar to the Fermi 
value at 100 MeV (Abdo et al. 2009b). Inserting values for L 7 from Table [21 one finds 
that star-forming and starburst galaxies could make a significant, > 10% contribution to the 
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EGB at 100 MeV, as previously suggested (Pavlidou & Fields 2002; Thompson et al. 2007). 

Observations with the Fermi Gamma-ray Space Telescope provide evidence that GeV 
emission has been detected from the starburst galaxy M82, and weaker though still significant 
evidence for detection of NGC 253. The Fermi LAT detections of these galaxies at GeV 
energies, together with the recent discovery of > 700 GeV 7 rays from M82 with VERITAS 
(Acciari et al. 2009) and > 220 GeV 7 rays from NGC 253 with H.E.S.S. (Acero et al. 2009), 
introduce a new class of 7-ray sources to 7-ray astronomy. Unlike 7-ray emitting blazars and 
radio galaxies powered by supermassive black holes, the evidence presented here supports 
a cosmic-ray origin for 7-ray production in starburst galaxies. Fermi observations over the 
upcoming years will improve our knowledge of spectra, variability properties, and number of 
7-ray bright starburst galaxies, which will also constitute important targets for observations 
with planned large Cherenkov telescope observatories CTA and AGIsB 
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Fig. 1. — Test statistic maps obtained from photons above 200 MeV showing the celestial 
regions (6° by 6°) around M82 and NGC 253. Aside from the source associated with each 
galaxy, all other Fermi-detected sources within a 10° radius of the best-fit position have been 
included in the background model as well as components describing the diffuse Galactic and 
isotropic 7-ray emissions. Black triangles denote the positions of M82 and NGC 253 at 
optical wavelengths; gray lines indicate the 0.68, 0.95, and 0.99 confidence level contours 
on the position of the observed 7-ray excess; green squares show the positions of individual 
background sources. The color scale indicates the point-source test statistic value at each 
location on the sky, proportional to the logarithm of the likelihood ratio between a 7-ray 
point-source hypothesis (Li) versus the null hypothesis of pure background (Lq); TS = 
2(lnLi - InLo) (Mattox et al. 1996). 



Table 1: Results of maximum likelihood analyses (gtlike) of M82 and NGC 253. 





RA a 


Dec a 


r 95 


F(> 100 MeV) b 


photon index 6 


significance c 




(deg) 


(deg) 


(deg) 


(10" 8 ph cm -2 s _1 ) 






M82 


149.06 


69.64 


0.11 


1.6±0.5 sta t ± 0.3 sys 


2.2±0.2 stat ± 0.05 sys 


6.8 


NGC 253 


11.79 


-25.21 


0.14 


0.6±0.4 stat ± 0.4 sys 


1.95±0.4 stat ± 0.05 sys 


4.8 



"Source localization results ( J2000) with r 95 corresponding to the 95% confidence error radius 
around the best-fit position. 

^Parameters of power-law spectral models fitted to the data: integrated photon flux > 100 

MeV and photon index. 

c Detection significance of each source. 
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Table 2: Properties of 7-ray galaxies lacking active central nuclei. 



Galaxy 


d 




M Gas 


F 7 a 


47rd 2 F 7 Q 






(Mpc) 


(yr- 1 ) 


(10 s Mq) 


(10~ 8 ph cm -2 s _1 ) 


(10 42 ph s" 1 ) 


(10 39 erg s" 1 ) 


LMC 6 


0.049 ± 0.001 


0.005 ± 0.002 


0.67 ± 0.08 


26.3 ± 4.7 


0.074 ± 0.013 


0.041 ± 0.007 


Milky Way c 


1 


0.02 ± 0.01 


6.5 ± 2.0 


4.6 ± 2.3 


5.5 ± 2.8 


3.2 ± 1.6 


M82 


3.6 ± 0.3 


0.2 ± 0.1 


2.5 ± 0.7 


1.6 ± 0.5 


25 ± 9 


13 ± 5.0 


NGC 253 


3.9 ± 0.4 


0.2 ± 0.1 


2.5 ± 0.6 


0.6 ± 0.4 


11 ± 7 


7.2 ± 4.7 



a 7-ray fluxes, F 7 , and luminosities, L 7 , computed in the energy range 100 MeV to 5 GeV. 

'LMC: distance measurement by Pietrzynski et al. 2009; supernova rate estimated by Tammann, et al. 1994; mass estimate by 
Bruns et al. 2005 (see also Westerlund 1997); 7-ray flux from Abdo et al. 2009a. 

c Gas mass estimate from Dame 1992; 7-ray flux from the Milky Way as viewed from a distance of 1 Mpc; 7-ray luminosity 
estimated using models which take into account pion-decay, inverse Compton, and bremsstrahlung photons produced in both 
the Galactic disk and halo (Bloemen et al. 1984, Strong et al. 2000, Pavlidou & Fields 2002). 




Fig. 2. — Relationship between supernova rate, total gas mass, and total 7-ray luminosity 
of four galaxies detected by their diffuse high-enery emission. In order of ascending 7-ray 
luminosity, the plotted galaxies are the LMC, Milky Way, NGC 253, and M82. Three panels 
are shown to compare different possible correlations with the 7-ray luminosity: total gas mass 
(left), supernova rate (center), and product of the total gas mass and supernova rate (right). 
This figure is based upon the observed quantities and associated uncertainties presented in 
Table [2 
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Fig. 3. — Spectral energy distributions of M82 and NGC 253. The spectra were obtained 
using gtlike, with flux points extracted based upon the parameters presented in Table [TJ 
Upper limits from the LAT correspond to the 0.68 confidence level. Three flux points in 
the TeV energy range are provided by VERITAS observations of M82 (Acciari et al. 2009). 
The single very high energy flux point for NGC 253 is computed from the integral photon 
flux over 220 GeV reported by the H.E.S.S. collaboration (Acero et al. 2009) and assumes a 
power-law spectral model marginalized over photon indices ranging from 2.0 to 3.0. Several 
theoretical predictions are plotted for comparison to the observed 7-ray spectra. 



